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Highly size-asymmetrical fluid mixtures arise in a variety of physical contexts, notably in suspensions of 
colloidal particles to which much smaller particles have been added in the form of polymers or nanoparticles. 
Conventional schemes for simulating models of such systems are hamstrung by the difficulty of relaxing 
the large species in the presence of the small one. Here we describe how the rejection- free geometrical 
cluster algorithm (GCA) of Liu and Luijten [Phys. Rev. Lett. 92, 035504 (2004)] can be embedded within 
a restricted Gibbs ensemble to facilitate efficient and accurate studies of fluid phase behavior of highly size- 
asymmetrical mixtures. After providing a detailed description of the algorithm, we summarize the bespoke 
analysis techniques of Ashton et al. [J. Chem. Phys. 132, 074111 (2010)] that permit accurate estimates of 
coexisting densities and critical-point parameters. We apply our methods to study the liquid-vapor phase 
diagram of a particular mixture of Lennard- Jones particles having a 10:1 size ratio. As the reservoir volume 
fraction of small particles is increased in the range 0-5%, the critical temperature decreases by approximately 
50%, while the critical density drops by some 30%. These trends imply that in our system, adding small 
particles decreases the net attraction between large particles, a situation that contrasts with hard-sphere 
mixtures where an attractive depletion force occurs. 



I. INTRODUCTION AND BACKGROUND 



Colloidal suspensions are a class of complex fluids that 
comprises systems as diverse as protein solutions, liq- 
uid crystals, and blood. Technologically, colloidal sus- 
pensions feature in applications such as coatings, pre- 
cursors to advanced materials, and drug carriersP One 
of the key issues in all these systems is the phase be- 
havior of the suspension, or more generally its stability. 
Attractive dispersion forces exist between uncharged col- 
loids that can engender phase separation or irreversible 
aggregation resulting in a gel — undesirable features in 
many applications. Accordingly, one seeks to control the 
phase behavior (as well as dynamical properties such as 
the rheology-2) by modifying the form of the effective 
interactions between the colloidal particles. There are 
several routes to achieving this, including charge stabi- 
lization (via modification of the pH) and steric stabi- 
lization (via grafting of flexible polymers onto the col- 
loidal surface) pEI Alternatively, the effective interactions, 
and hence colloidal phase behavior, may be manipulated 
through the addition of nanoparticles, with nanoparticle 
size, concentration, and charge as control parametersPHS 
The simplest and most celebrated example concerns col- 
loids which interact (to a good approximation) as hard 
spheres. Adding nanoparticles in the form of small non- 
adsorbing polymers engenders an attractive "depletion" 
force between the colloidal particles.— This attraction can 
drive phase separation resulting in a colloid-rich ('liquid') 
phase and a colloidal-poor ('gas') phase^ — a phenomenon 
akin to the fluid-fluid transitions occurring in molecular 
liquids and their mixtures. Yet richer behavior occurs 
when one transcends simple hard-sphere potentials be- 



tween the nanoparticles and the colloids. For example, if 
the nanoparticles are weakly attracted to the colloids but 
repel one another, they can form a diffuse (nonadsorbed) 
"halo" around each colloid particle! 5 ! 9 " 13 ! The net influ- 
ence on the effective colloid-colloid interaction depend s 
on the nanoparticle density in a nontrivial waypEflHI 

In view of the broad range of effects that can arise 
when nanoparticles are added to a colloidal suspension, 
their prototypical model representation, namely a size- 
asymmetric fluid mixture, has attracted considerable the- 
oretical and computational attention over the past years. 
Analytical approaches typ ically either focus on drasti- 
cally simplified models-^ or attempt to render the size 
asymmetry tractable by integrating out the degrees of 
freedom associated with the small (nano)particles (see 
Ref. [TS] for a review) . The latter strategy yields a one- 
component system of colloids described by an effective 
pair potential representing the net influence of the small 
particles. One shortcoming of this approach is that, for 
all but the simplest types of nanoparticles, the map- 
ping to a one-component system is approximate, be- 
cause it neglects many-body colloidal interactions that 
can considerably alter the nanoparticle distributions and 
hence the interactions induced by them. These effects 
may be significant in the regimes of density at which 
phase separation occurs.^ Recent work has addition- 
ally raised concerns rega rding the accuracy of effective 
potentials in this regimepQH and has also emphasized 
the significance of corr ection s to the entropic depletion 
picture in real colloid s^ 9 -* 2 ^ as well as the importance 
of polydispersityp322l On the other hand, various com- 
putational techniques, most notably Monte Carlo (MC) 
methods, are capable of explicitly incorporating fluctua- 
tion and correlation effects. Conventional MC techniques 



(which attempt to displace or insert and delete particles) 
are restricted to fluids mixtures in which the size ratio 
is of order unity, rendering them unsuitable for the sim- 
ulation of colloid-nanoparticle suspensions, where typ- 
ical size ratios encountered can extend to one or two 
orders of magnitude.^ The computational bottleneck re- 
sults from the effective "jamming" of the large species by 
even low volume fractions of the small particles. How- 
ever, this problem has been resolved by means of th e ge- 
ometric cluster algorithm (GCA) of Liu and LuijtenpSES 
in which configuration space is sampled via rejection-free 
collective particle updates, each of which facilitates the 
large-scale movement of a substantial subgroup of parti- 
cles (a "cluster"). Although the original algorithm op- 
erates in the canonical ensemble, and hence cannot ad- 
dressphase separation phenomena directly, in previous 
worhpS W e have developed a generalization that embeds 
the GCA in the restricted Gibbs ensemble (RGE), such 
that clusters containing both large and small particles are 
exchanged between two simulation boxes of fixed equal 
volumes. The resulting density fluctuations within one 
box can be analyzed to determine the phase behavior. 

The purpose of the present paper is firstly to provide 
a more detailed description of the basic GCA-RGE al- 
gorithm that was introduced in Ref. 1251 and to integrate 
recent advances that we have made in data analysis meth- 
ods for determining coexistence and critical-point prop- 
erties within the RGEP3 We then apply the improved 
methodology to study the liquid-vapor coexistence prop- 
erties of a mixture of Lennard- Jones (LJ) particles having 
a size ratio q = 0.1. In so doing we adopt one aspect of 
effective fluid approaches, namely we focus on the liquid- 
gas phase coexistence properties of the large species (col- 
loids) which are assumed to be immersed in a supercriti- 
cal fluid of small particles of quasi-homogeneous density. 
This choice of perspective mirrors the experimental real- 
ity, namely that often only the colloidal particles can be 
individually imaged. Accordingly, the phase diagrams 
that we present are single-component projections (i.e., 
referring to the large species) of the full phase diagram, 
obtained at a prescribed reservoir volume fraction of the 
small species, which we vary in the range 0-5%. Note, 
however, that the small particles are treated explicitly 
and exactly in our simulations, making this method su- 
perior to effective-potential approaches which integrate 
out the degrees of freedom associated with the small par- 
ticles. 

This paper is arranged as follows. Section |ll] intro- 
duces our model system, a binary Lennard- Jones fluid. 
The GCA-RGE MC algorithm capable of simulating this 
system in the highly size-asymmetrical limit is described 
in Sec. |IH| together with an outline of techniques for de- 
termining phase coexistence properties and critical-point 
parameters within the RGE. Moving on to our results, 



Sec. IV presents measurements of the large-particle co- 
existence densities as a function of the reservoir volume 
fraction of small particles. We also discuss the under- 
lying reasons for the observed trends in the coexistence 



properties in terms of measurements of the fluid struc- 
ture. Finally, Sec. fV] considers the implications of our 
findings, the efficiency of our simulation approach com- 
pared to more traditional schemes, and an outlook for 
further work. 



II. MODEL SYSTEM 

The model with which we shall be concerned is a bi- 
nary mixture of spherical particles, whose two species are 
denoted I (large) and s (small) . Pairs of particles labeled 
i and j (having respective species labels ji and jj) inter- 
act via a Lennard- Jones potential, 



<M0=4e 



iiij 



( u ini \ _ ( u ini \ 



(1) 



where e 7i7 is the well depth of the interaction and 



<7 7i7j sets its range based on the additive mixing rule 
<7 7i7j = (<r 7i +<r 7i )/2. o~ 1i and a lj represent the particle 
diameters. Interactions arc truncated at r 
and we take o\ as our unit length scale. 
In Sec. IV we study the case q = a ss /au 
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10:1 size ratio. We shall determine the phase coexistence 
properties of the large particles as a function of tempera- 
ture for a prescribed reservoir volume fraction 77J of small 
particles, as controlled by the imposed chemical potential 
of small particles fj, s . It should be noted, however, that 
since the small particles are not infinitely repulsive, their 
volume fraction is notional in the sense that we use the 
value of a s as if it were a hard-core radius, i.e., we take 
rf s = TrN s a^/(6V) where N s is the average of the fluc- 
tuating number of small particles contained within the 
system volume V. 

Since we adopt the viewpoint that the small par- 
ticles act as a background to the large ones, we set 
e ss = eis — e«/10, which ensures that the small-particle 
reservoir fluid is supercritical in the temperature range 
of interest here, namely down to well below the critical 
point of the large particles. It is therefore natural to de- 
fine the dimensionless temperature T in terms of the well 
depth of the interaction between the large particles, i.e., 
T = k^T/su, where fee is Boltzmann's constant and T 
the absolute temperature. 



III. METHODOLOGY 

A. The GCA-RGE algorithm 

In the original GCAp2 a fixed number of particles is 
located in a single, periodically replicated simulation box 
of volume V. These particles are then moved around via 
cluster moves, in which a subset of the particles (iden- 
tified by means of a probabilistic criterion) is displaced 
via a geometric symmetry operation. To realize density 
fluctuations, we employ two simulation boxes and ex- 
change particles between both boxes, as in the Gibbs 



ensembleP^ However, rather than exchanging individual 
particles, we use the GCA to exchange entire clusters of 
particles, so that we retain the primary advantage of the 
GCA, namely the rapid decorrelation of size-asymmetric 
mixtures. As in the original GCA, a variety of sym- 
metry operations is possible; to connect to the original 
description,^ we phrase the algorithm here in terms of 
point reflections with respect to a pivot. Since a point 
reflection will generally displace some particles outside 
of the original simulation cell, we need to adopt periodic 
boundary conditions for both simulation cells. Moreover, 
as will transpire below, all particles that belong to a clus- 
ter and that are part of the same simulation cell will 
retain their relative positions during the cluster move. 
Thus, the two simulation cells must have the same di- 
mensions. This symmetric choice, in which both cells 
have an identical, constant volume V, is referred to as 
the RGE. 

We first describe the GCA-RGE for the case of a single 
species of particles that interact through an isotropic pair 
potential V(r). N = Ai + N 2 particles are distributed 
over the two simulation cells, with N\ particles in simula- 
tion cell 1 and A2 particles in simulation cell 2. Nq is cho- 
sen to match a desired average density po — N /(2V). A 
cluster move within the GCA-RGE proceeds as follows. 
A pivot is chosen at a random position within simulation 
cell 1 and a second pivot is placed at the corresponding 
position within simulation cell 2. One of the No parti- 
cles is chosen as the seed particle of the cluster. This 
particle i, which thus can be located in either simulation 
cell, is point-reflected with respect to the pivot (in its 
own simulation cell) from its original position r^ to the 
new position r^. However, rather than placing the par- 
ticle at the new position (modulo the periodic boundary 
conditions) in its original box, we place it at the corre- 
sponding position y[ in the other box. Subsequently, in 
keeping with the methodology of the GCA, all particles 
in the first box that interact with particle i in its original 
position (the "departure site") r^ as well as all parti- 
cles in the second box that interact with particle i in its 
new position (the "destination site" ) r \ are considered for 
point reflection with respect to the pivot point in their 
respective box and subsequent transfer to the opposite 
box. These particles, which we refer to with the index j, 
are point-reflected and transferred with probability 



Pij = max[l - exp(-/3Ajj), 0] 



(2) 



where = l/(k&T) and Ay = — V(|r^ — r^ |) if i and j 
reside (prior to the transfer of particle i) in the same 
cell. If i and j initially reside in different cells (and 
hence do not interact prior to the transfer of particle i), 
Ajj = V(\f' t — Tj\). This process is repeated iteratively, 
i.e., for each particle j that is transferred to the oppo- 
site box, all neighbors that interact with j either near its 
departure site or near its destination site, and that have 
not yet been transferred in the present cluster step, are 
considered for point reflection and transfer as well. This 
process proceeds until there are no more particles to be 



considered; all particles that are indeed point-reflected 
and transferred are collectively referred to as the cluster. 
Observe that the pair energy of all particles that are part 
of the cluster remains unchanged: If two particles reside 
in the same simulation cell prior to the cluster construc- 
tion and both become part of the cluster, their separation 
remains constant. Likewise, if two particles reside in dif- 
ferent cells prior to the cluster construction and both are 
transferred, then their interaction energy is zero before 
and after the cluster move. The same holds true for the 
pair interactions between all particles that are not part 
of the cluster. Thus, the total energy change induced by 
the cluster move originates from the change in pairwise 
interactions between members of the cluster and parti- 
cles that are not part of the cluster. In the terminology 
of Ref. 23 such "bonds" are either broken if a particle is 
included in the cluster whereas a neighbor near its depar- 
ture site is not, or formed if a particle interacts with a 
neighbor near its destination site, and this neighbor does 
not become part of the cluster. 

Although the cluster formation process is probabilis- 
tic, we note that pij only depends on the pair poten- 
tial between particles i and j, rather than on the to- 
tal energy change resulting from the displacement and 
transfer of particle j. As a result, the cluster algo- 
rithm is self-tuning: Overlaps of repulsive particles will 
be avoided and strongly bound particles tend to stay to- 
gether. Indeed, owing to the choice of the bond probabil- 
ity Pij: Eq. (|2|), no further acceptance criterion needs to 
be applied upon completion of the cluster, leading to a 
rejection-free algorithm in which large numbers of parti- 
cles are moved nonlocally. The proof of detailed balance 
is identical to that provided in Ref. 24 for the original 
GCA, where it was demonstrated that the ratio of the 
probability of constructing a cluster in a given config- 
uration X leading to a configuration Y [the transition 
probability T(X -> Y)} and the reverse transition proba- 
bility T(Y —> X) is the inverse of the ratio of Boltzmann 
factors of the respective configurations. The presence of 
two simulation cells simplifies rather than complicates 
the proof, just like Ajj in Eq. pi is a special case of the 
original expressionpS owing to the fact that two particles 
do not interact if they reside in opposing boxes. 

The generalization to multiple species is straightfor- 
ward, and does not lead to any conceptual changes in 
the algorithm. Indeed, the GCA shows its primary ad- 
vantages in the simulation of size-asymmetric mixtures, 
as it realizes nonlocal moves without the usual decrease in 
acceptance ratioPS However, whereas there is no limita- 
tion on the number of species, the overall volume fraction 
must be kept below a threshold value. Above this thresh- 
old, which is related to the percolation threshold and de- 
pends on system composition and interaction strengths 
between the particlesp*! the cluster frequently contains 
the majority of all particles. This is detrimental to the 
performance of the algorithm, as it is computationally 
expensive to construct such clusters, whereas the config- 
urational change in the system is very small. The exis- 



tence of this threshold also necessitates the use of an im- 
plicit solvent, as is common in the simulation of colloidal 
suspensions. Moreover, for reasons explained in detail in 
Sec. |III C[ in our simulations of binary mixtures we com- 
bine the geometric cluster moves with grand-canonical 
moves for the small species. It is important to empha- 
size that the different types of MC moves are indepen- 
dent. Thus, the small species fully participate in the 
cluster construction process and the advantage of non- 
local reject ion- free moves is retained, yet the density of 
small particles in both simulation cells is controlled by a 
chemical potential p s . 

In Ref. 1241 a number of technical improvements to 
the GCA are described. These can all be applied to 
the GCA-RGE algorithm. Most notably, it is possible 
to decrease the average cluster size, and hence increase 
the packing fractions that can be simulated efficiently, 
through biased placement of the pivot. Furthermore, for 
mixtures of particles with large size disparities, the clus- 
ter construction process can be facilitated by employing 
multiple subcell structures and corresponding neighbor 
Ms.™ 

Lastly, we note that, during the preparation of our 
original workp*! BuholP^ proposed an approach that has 
significant similarities to the GCA-RGE method. His 
method also employs two boxes of identical size, and ex- 
changes clusters of particles. However, rather than the 
GC Apg kg ugeg ^Yxq original geometric algorithm of Dress 
and KrauthpSl which is only applicable to hard spheres. 
Moreover, for each point reflection it is decided at random 
whether a particle is transferred to the opposing box or 
not. If this decision were made only once per cluster (i.e., 
upon selecting the seed particle of the cluster) , this would 
amount to an alternation of the GCA-RGE with regular 
GCA moves. On the other hand, if it is decided indepen- 
dently for each particle that is added to the cluster, the 
average cluster size will be larger than in the GCA-RGE, 
generally an undesirable situation. The most important 
difference, however, between, our approach and Ref. HH1 
is that the latter can only be used for the idealized case 
of symmetric binary mixtures, where the critical compo- 
sition is known a priori. By contrast, in our method we 
employ the relationship between the RGE and the grand- 
canonical ensemble to derive a prescription for locating 
the critical point and coexistence curve for general binary 
mixtures. 



B. Locating phase coexistence and criticality in the RGE 

The absence of volume exchanges between both simula- 
tion boxes in the symmetrical restricted Gibbs ensemble 
implies that, unlike for the full Gibbs ensemblep-^ there 
is no automatic pressure equality and hence no guarantee 
that the measured particles densities are representative of 
coexistence. In this section we outline how one can nev- 
ertheless extract coexistence properties from RGE simu- 
lations without resorting to direct measurements of pres- 



sure. A fuller account of the theoretical basis of the meth- 
ods we describe can be found in Ref. 26. 

Within the RGE framework for our mixture, the to- 
tal density po of large particles across the two boxes is 
fixed. However, the one-box density of large particles, 
p = Ni/V, fluctuates. For any given choice of po, the 
form of the probability distribution of p, Pi, (p) , depends 
both on the temperature T and on the choice of the 
chemical potential p s of the small particles. As shown 
in Ref. 26, measurement of the form of Pl(p) for a range 
of values of po provides a route to the coexistence and 
critical-point parameters. The basic strategy is as fol- 
lows. Within the RGE, one explores the coexistence 
region by varying po at fixed T and p s . For po suffi- 
ciently far inside the coexistence region, the distribution 
Pl(p) exhibits a double-peaked form, with peaks located 
at densities p_ and p + . In general, however, these peak 
densities do not coincide with the gas and liquid coex- 
istence densities p gas and pu q — a situation which con- 
trasts with the full Gibbs ensemble. An important excep- 
tion is when po equals the coexistence diameter density 
p d = (p gas + pu q )/2, for which one findsPsl 



P- = Pgas 
P+ = Pliq 



when po = p d 



(3) 



Another important case is when po = (pgas + Pd)/2 for 
which one finds 



P- = Pgas 

p+ = Pd 



when po = {p gas + Pd)/2 . (4) 



Enforcing consistency between Eqs. (pi) and Q suffices 
to permit determination of p c i and hence [via Eq. p|] 
the coexistence densities. It is convenient to achieve this 
graphically (see Fig. [I]) by plotting the low density peak 
p_ both against p and against 2p — p_ : The value of po 
at which the two curves intersect provides an estimate for 
the coexistence diameter pd and one can simply read off 
the coexistence densities from the corresponding values 
of p_ and p + . In Ref. |2"61 this "intersection method" was 
shown to be very accurate for determining coexistence 
properties and to exhibit finite-size effects comparable 
to those found in grand-canonical simulations. Indeed it 
turns out to be much more accurate than the technique 
we proposed previously for determining the coexistence 
diameter in the RGE, 25 wherein one determines pd as the 
value of po at which the variance of Pl{p) is maximized. 
We have found this latter procedure to be considerably 
more sensitive to finite-size effects than the intersection 
method, and it was therefore not used here. 

Turning now to the matter of estimating critical pa- 
rameters within the RGE, an accurate technique for 
achieving this has been described in detail in Ref. [251 
The basic idea is to measure the "iso-Q* curve" intro- 
duced in Ref. [25j which is simply the locus of points in 
Po~T space for which the fourth-order cumulant ratio of 
Pl(p), 
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FIG. 1. (a) Illustration of the operation of the intersection 
method described in the text for the determination of the co- 
existence diameter density. Data is shown for the state point 
rf s = 0.01, T = 0.8(= 0.764T C ). Plotted are measured es- 
timates of the average value p_ of the low-density peak of 
Pl(p\po), for a series of values of po. The same data is also 
shown plotted against 2po — p- ■ The value of po at which the 
two data sets intersect [po = 0.351(3)] serves as an estimate 
of the coexistence diameter density pd- (b) and (c): The mea- 
sured peaks of Pl(p) for po = pd, whose individual integrated 
averages yield estimates of the coexistence densities. 



__ ((P-P0) 2 ) 2 
" ((P-Po) 4 ) 



(5) 



matches the independently knowrPSl fixed-point value 
Q* = 0.711901 appropriate to the Ising universality class 
and the RGE ense mbleP^ 

Now, it transpire d 25 * 26 ! that the iso-Q* curve is essen- 
tially a parabola in p$-T space, the position of whose 
maximum represents a finite-size estimator of the critical- 
point density p c and temperature T c . This maximum can 
be accurately located via a simple quadratic fit to mea- 
sured points on the iso-Q* curve. In practice we deter- 
mine this curve as follows. A simulation is performed 



FIG. 2. Estimates of points on the iso-Q* curve for rf s — 0.01, 
obtained for a system of size L = 10, as described in the 
text. A parabolic fit to the data (solid line) identifies the 
coordinates of the maximum of the curve which serves as a 
finite-size estimate of the critical-point parameters. 



at some po and T to determine P{p). This distribu- 
tion is then extrapolated in temperature via histogram 
reweighting22 to determine that temperature for which 
Q = Q* . The procedure is then repeated for a range of 
values of p allowing us to trace out the whole iso-Q* 
curve. An example of the resulting form of this curve 
is shown in Fig. [2j Note that, in general, for reasons of 
computational economy, the majority of the points that 
we determine on an iso-Q* curve are for densities lower 
than the critical density, since the efficiency of the cluster 
algorithm is greater at lower overall volume fractions of 
particles. 

Estimates of the critical parameters obtained from the 
iso-Q* maxima for a range of system sizes can, in prin- 
ciple, be extrapolated to the thermodynamic limit using 
finite-size scaling relations derived in Ref . [55] which fully 
account for both field-mixing effects and corrections to 
scaling. Unfortunately, in the present work, the com- 
putational cost of simulating more than one system size 
was found to be prohibitive. However, the variations that 
we find in critical-point parameters as a function of rf s 
dwarf those that one might expect on the basis of finite- 
size effects alone. Thus we are nevertheless able to report 
reliable trends from our measurements. 



C. Treatment of the small particles 

As was argued in Sec. [TJ it is the relaxation of the 
large particles that constitutes the sampling bottleneck 
for highly size-asymmetrical mixtures. Local MC up- 
dates of small particles are computationally relatively 
unproblematic, irrespective of whether one performs par- 
ticle displacements or insertions and deletions. Conse- 
quently, one has the choice of treating the small particles 



canonically so that their density is globally conserved, 
or grand canonically, in which case the density fluctuates 
under the control of a prescribed chemical potential. The 
choice one makes in this regard greatly affects the man- 
ner in which the bulk phase behavior is probed. Moreover 
it transpires that only the grand-canonical treatment of 
the small particles is compatible with our intersection 



method (Sec. Ill B ) for determining coexistence parame- 
ters. 

To clarify these points, we show in Fig. [3] sketches of 
the isothermal bulk phase diagram of an exemplary size- 
asymmetrical binary mixture with large-particle den- 
sity pi, small-particle density p s , and conjugate chemi- 
cal potentials pi and p s , respectively. Figure p[ai) shows 
the phase behavior in the pi~p s plane, while the corre- 
sponding phase diagram in the pi~p s plane is shown in 
Fig. pJtaii). In constructing these sketches we have antic- 
ipated the behavior of the model of Sec. [TTJ namely that 
the larger species has stronger attractive interactions and 
thus phase separates on its own [vertical axis of Fig. [3l ai)] 
at the chosen temperature, while the small-particle fluid 
(horizontal axis) does not. With interaction strengths 
chosen in this way, the larger particles will typically ac- 
cumulate in the liquid phase, with its shorter interparticle 
distances, as shown by the representative tie lines in the 
density representation's Note that in chemical-potential 
space, coexistence occurs on a line of points, as shown in 
Fig.gaii). 

Let us consider first the canonical scenario in which one 
traverses the coexistence region at constant bulk density 
of small particles, as expressed by the dashed trajectory 
included in Figs. [3lai) and (aii). Clearly the tie lines 
in Fig. p[ai) cross any line of constant p s moving from 
smaller values of pi/p s at the gas end to larger ones for 
the coexisting liquid. Hence this path generates a se- 
quence of pairs of coexistence states, one for each tie line 
crossed. The same trajectory in terms of the chemical 
potentials pi~p s is shown in Fig. kMaii). Here the path 
followed first meets the coexistence line, tracks along it 
it for some distance and then separates from it. 

The alternative (grand-canonical) scenario, in which 
the small-particle density is permitted to fluctuate at 
constant chemical potential p s , is illustrated in Figs. 
p[bi) and (bii). Here, as pi is varied at constant p s , 
the system crosses the coexistence line at a single point. 
In density space the corresponding trajectory thus fol- 
lows a particular tie line though the coexistence region, 
as shown in Fig. [3lbi) . 

In seeking to apply the RGE ensemble to study a bi- 
nary mixture, it is therefore imperative that one adopts a 
grand-canonical treatment of the small particles. Doing 
so ensures that only a single pair of coexistence states is 
encountered inside the bulk coexistence region, i.e., that 
one tracks a tie line of the bulk phase diagram. This is a 
prerequisite for the correct operation of our intersection 
method, which is designed to determine first the diameter 
density for a single pair of coexistence states as a prelude 
to determining the coexistence densities of the large par- 
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FIG. 3. Isothermal cuts through the exemplary phase diagram 
of a binary fluid mixture described in the text. Liquid-vapor 
coexistence is represented in terms of (i) densities (pi-p B ) and 
(ii) chemical potentials (pi-p a ). In (a) the coexistence region 
is crossed along a path of constant p s , whereas in (b) it is 
crossed along a path of constant p a - 



tides themselves. We further note that a grand-canonical 
treatment of the small particles corresponds more closely 
to common experimental arrangements where one typi- 
cally measures properties of the mixture with respect to 
variations of a reservoir volume fraction of small parti- 
cles. 



IV. RESULTS 

Equipped with the methods described above, we can 
set about the task of determining the coexistence proper- 
ties of the large particles in the presence of a sea of small 
ones. To this end we apply the GCA-RGE method to 
study liquid- vapor phase coexistence in a q = 0.1 LJ mix- 
ture (see Sec. III]). In addition to the cluster updates which 
swap whole groups of particles (including both large and 
small species) between boxes, small particles are sampled 
across both boxes using a standard local grand-ca nonica l 
algorithm at constant chemical potential (see Sec. |III C ) . 
As discussed above, we choose p s to yield (for each tem- 
perature of interest) a prescribed volume fraction if s of 
small particles in the reservoir. This requires prior knowl- 
edge of the reservoir equation of state rf s {p s , T), which we 
obtained via explicit simulation of the pure fluid of small 
particles. Note that the computational cost of obtaining 
rf s (p Sl T) is low, particularly if one employs histogram 
extrapolation^ to scan a region of p and T surrounding 
each simulation state point. 

The GCA-RGE simulations are performed using two 
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FIG. 4. Distribution of the fraction of large particles in the 
cluster, / C) at criticality for the values of rf s shown in the 
legend. 



cubic periodic simulation boxes of linear size L = 10. 
We consider seven values of the reservoir volume fraction 
of the small particles, rf 3 = 0.005, 0.01, 0.015, 0.02, 0.03, 
0.04, and 0.05. In the limit of low densities of large parti- 
cles, these values of r] r 8 correspond to average numbers of 
small particles in the range 10 4 -10 5 . The computational 
expenditure incurred in simulating such great numbers 
of small particles places an upper bound on the value of 
rf s for which it is feasible to perform a full determina- 
tion of the coexistence binodal. Further difficulties arise 
from the fact that the typical cluster size at coexistence 
was found to grow steadily as we increased rf s . This is 
demonstrated in Fig. [4] which plots the distribution of the 
fraction of large particles in the cluster for rf s = 0.01, 0.03 
and 0.05 measured at the respective critical point param- 
eters. These distributions are bimodal, with some fairly 
small clusters comprising just a few particles, and many 
clusters that comprise the vast proportion of large par- 
ticles. Updating such clusters results in only relatively 
minor alterations to a configuration and consequently, 
we are able to determine the coexistence binodal only for 
r/s < 3%, whereas for r]l = 0.04 and 0.05 we restrict our- 
selves to determining critical-point parameters. It should 
be stressed however, that the cluster sizes observed in 
the present study may not provide a general guide to the 
maximum rf s at which the GCA-RGE scheme will oper- 
ate. This is because, as we shall show, our choice of in- 
terspecies interactions engenders a large depression in T c 
with increasing nf s which in turn promotes the formation 
of large clusters due to the temperature dependence of 
the GCA bond- formation probability Eq. (pi). However, 
other choices of interactions can be expected to lead to 
a different temperature dependence of the critical point 
parameters, hence allowing larger values of rf s to be at- 
tained. 

The critical-point parameters are determined, for each 
rf s studied, from measurements of the iso-Q* curve. For 



FIG. 5. Phase diagrams showing the liquid (diamonds) and 
gas (squares) coexistence densities of large particles for q = 
0.1. Data are shown for reservoir volume fractions (top to 
bottom) rf s = 0, 0.005, 0.01, 0.02, and 0.03. Also shown 
in each case is the coexistence diameter (circles) and critical 
point (asterisks). 
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rf s < 3%, the intersection method described in Sec 
is deployed to determine the large-particle coexistence 
densities in the subcritical regime. Figure [5] presents our 
results for the p-T binodal. The principal feature is - 
as previously mentioned - a strong depression of the bin- 
odal to lower temperatures and lower densities as rf s is 
increased. The scale of the associated shifts in the crit- 
ical parameters is made apparent in Fig. [6j which plots 
our estimates of the critical temperature and density as a 
function of rf s . One sees that as the reservoir volume frac- 
tion of small particles is increased in the range - 0.05, 
the critical temperature decreases by approximately 50%, 
while the critical density drops by some 30%. The error 
bars shown on the estimates of the critical parameters 
derive from a bootstrap analysis of the various quadratic 
fits that are consistent with the uncertainties in the locus 
of each iso-Q* curve (Fig. [7]). 

To demonstrate the correctness of our method, we 
compare the binodal for rf s — 0.01 with that obtained 
using a quite different approach, recently proposed by 
two of us.^l This is a fully grand-canonical MC scheme 
in which large particles are gradually transferred to and 
from the system by means of staged insertions and dele- 
tions. To negate ensemble differences that occur when 
comparing results in finite-size systems, we transform 
the grand-canonical distribution of the large-particle den- 
sity, P(p), to the RGE using the exact transformation 
P(p) = P(p)P(2p Q - p ) .1251261 We tnen proceed to locate 
coexistence as if the data had been generated in the RGE 
by treating po as a parameter of the transformation. The 
resulting coexistence densities are compared with those 
obtained via the GCA-RGE simulations in Fig. [Hj The 
agreement is good, particularly at low temperature. The 
deviations near criticality arise from the difference in the 
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FIG. 6. (a) Critical temperature and (b) critical density ver- 
sus rf s as determined from the iso-Q* curves. Error bars derive 
from a bootstrap analysis with 100 resamples. 
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FIG. 7. The measured iso-Q* curves for (top to bottom) 
rf s = 0, 0.005, 0.01, 0.015, 0.02, 0.03, 0.04, and 0.05. Also 
shown are the estimated error bars from which we calculated 
the overall uncertainty in critical parameters via a bootstrap 
analysis. 



FIG. 8. Comparison for rff. — 0.01 of the binodal obtained 
using the GCA-RGE technique (crosses) and the grand- 
canonical approach (circles) described in Ref. 1341 Statistical 
uncertainties are comparable to the symbol sizes. Differences 
in the results near criticality arise from the different system 
sizes used in the two cases. 



system size used in each case (L = 7.5au for the grand- 
canonical system and L = 10er;; for the RGE system), 
and thus reflect that the correlation length exceeds the 
system size in the grand-canonical simulation. 

It is instructive to attempt to relate the shift in the 
binodal occurring with increasing small-particle density 
to alterations in the underlying local fluid structure. An 
indication as to the factors at work here follows from 
a study of the effect of small particles on the effective 
potential between a pair of large particles, 



/3W(r) = lim - 

p->0 



Hmi( r )} 



(6) 



An example is shown in Fig.pMa) which compares f3W(r) 
for the cases rf s = (which simply corresponds to the 
bare Lennard- Jones potential) with that for rf s = 0.05, 
at T = 1.3. One sees that for typical separations of large 
particles, the effective potential is less attractive than the 
bare interaction. Thus the net effect of the small particles 
is repulsive as shown by the difference plot in fig. pMa), 
a feature that accords with the reduction in the critical 
temperature. A likely reason for this is to be found in 
the associated form of gi s (r) describing the correlations 
between a large particle and a small particle, as shown 
in Fig. |9jb) at r/l — 0.05. This shows that small particles 
form a diffuse, non absorbing cloud around each large 
particle because of their weak mutual attraction. Pre- 
sumably, however, the free-energy cost arising when the 
clouds associated with two or more large particles over- 
lap acts to reduce the intrinsic attractions between large 
particles. Interestingly, the difference plot of Fig. [9]ja) 
shows that at very small separations of large particles 
(corresponding to high overlap energy) the effect of the 
small particles changes from being repulsive to being at- 
tractive. 
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FIG. 9. (a) The measured form of the effective potential 
PW(r) defined in the text, at temperature T — 1.3. Data 
are shown for the bare LJ potential (rf g — 0, dashed line) 
and rf s = 0.05 (solid line) and their difference (dashed-dotted 
line), (b) Form of gi s {r) for p n ->• at 77J = 0.05, T = 1.3. 




FIG. 10. Configurational snapshot of the two boxes in the 
restricted Gibbs ensemble at coexistence for rf s — 3%, T = 
0.88T C . 



Finally, we show in Fig. [10] a configurational snapshot 
of our simulation boxes at coexistence (i.e., po = pd) for 
the case rf s = 3%, T — 0.88T C . This provides a visual im- 
pression of the character of the coexisting phases and the 
extent to which the large particles are severely "jammed" 
by the small ones. 



V. DISCUSSION AND CONCLUSIONS 

In summary, we have described a variant of the Geo- 
metric Cluster Algorithm 23 for the accurate determina- 
tion of phase behavior in highly size-asymmetrical fluid 
mixtures. The method (an early version of which was 
previously described in Ref. H5J operates by swapping 
clusters containing large and small particles between two 
boxes of equal volume, the global density of large par- 
ticles being fixed. The resulting spectrum of single-box 
fluctuations of the large particles can be analyzed with 



respect to changes in their global density using the in- 
tersection method of Ashton et alW* to yield accurate 
estimates of coexistence densities. Critical points can 
similarly be located to high precision by using an ap- 
propriate finite-size estimator for criticality, namely the 
maximum of the iso-Q* curve. 

We have applied the method to a LJ mixture with size 
ratio 10:1 to determine the coexistence properties of large 
particles for small-particle reservoir volume fractions in 
the range < rf s < 3%. Additionally, critical-point pa- 
rameters were determined for rf s = 0.04 and 0.05. Our 
results show that when the small particles are weakly 
attracted to the large ones, their net effect is to lower 
the degree of attraction between large particles. As a 
consequence, the coexistence binodal shifts to lower tem- 
peratures, confirming the preliminary findings of Ref. 1251 
Such a situation contrasts markedly with the depletion 
effect applicable to small particles that interact with the 
large ones like hard spheres,^ for which there is a net 
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increase in the degree of attraction between large parti- 
cles. Our measurements of local structure suggest that 
in the case we have considered, the small particles form a 
diffuse (nonadsorbing) cloud surrounding each large par- 
ticle. The overlap of clouds necessary for two large par- 
ticles to approach one another appears unfavorable in 
free-energy terms, leading to a net decrease in the degree 
of attraction between large particle s. T his is reminiscent 
of the "nanoparticle haloing" effect P^ 

It is gratifying to note that the GCA-RGE method 
permits the study of phase behavior in regimes that are 
inaccessible to traditional simulation approaches. Specif- 
ically, the phase diagrams we have presented could not 
have been obtained using even the most efficient tradi- 
tional approach to fluid phase equilibria, namely stan- 
dard grand-canonical simulation!^ For instance, for rf s > 
0.005 our tests show that the grand-canonical relaxation 
time is too large to be reliably estimated. Nevertheless, 
a lower bound on the grand-canonical relaxation time, 
relative to that of the pure LJ fluid, can be estimated 
via a comparison of the large-particle transfer (inser- 
tion/deletion) acceptance probability p acc . For liquid- 
like densities of the large particles (p « 0.6), p acc is of 
the order of 10 -4 at tf s = 0.005. 25 Upon increasing rjl 
to a volume fraction of merely 1%, this probability falls 
to p acc ~ 10 -6 . These values are to be compared with 
Pace ~ 10 _1 for the pure LJ fluid. One can therefore ex- 
pect the grand-canonical relaxation time of the mixtures 
studied here to be several orders of magnitude greater 
than for the pure LJ fluid. 

Notwithstanding the efficiency gains provided by the 
GCA-RGE approach, it should be stressed that the re- 
sults we have reported nevertheless entailed a significant 
computational outlay. Specifically, runs to determine 
each coexistence point typically varied in length between 
100 and 3,000 hours of CPU time on a 3 GHz proces- 
sor. The upper value in this range was that required at 
the highest volume fraction of small particles studied (for 
which there are very many small particles) and the low- 
est temperature (where most of the large particles are 
involved in each cluster update). Thus whilst studies of 
phase behaviour in highly asymmetrical mixtures cannot 
yet be regarded as routine, they are now at least feasible. 

With regard to future studies of highly asymmetrical 
mixtures, one barrier to attaining higher values of rf s 
and smaller values of q is simply the computational over- 
head associated with large numbers of small particles, 
although we note that the GCA has been successfully 
applied to systems with millions of nanoparticles.-^ Ad- 
ditionally, in the present model, the suppression of the 
critical temperature with increasing 77^ leads to a rapid 
growth in the cluster size which renders the GCA-RGE 
approach increasingly less efficient. More generally, how- 
ever, in situations where the small particles induce an 
effective (depletion) attraction between the large ones, 
we expect that the cluster size will remain manageable 
to rather larger r]l than studied here. 

Finally, we mention an alternative approach, proposed 



by two of us, for determining coexistence properties in 
highly size-asymmetrical mixtures!^ This method uti- 
lizes an expanded grand-canonical ensemble in which the 
insertion and deletion of large particles is accomplished 
gradually by traversing a series of states in which a large 
particle interacts only partially with the environment 
of small particles. Implementing this approach requires 
prior determination of a multicanonical weight function 
to bias insertions of the particles, and thus renders it 
less straightforward to use than the GCA-RGE. How- 
ever, being fully grand canonical does have the advan- 
tage of providing information on the chemical potentials 
of large particles, thereby permitting histogram reweight- 
ing in terms of density as well as temperature. In future 
work we hope to provide a systematic comparison of the 
relative computational cost of both approaches in various 
parameter regimes. 
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